Resonant-Cavity-Induced Phase Locking and Voltage Steps in a Josephson Array 
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We describe a simple dynamical model for an underdamped Josephson junction array coupled to a 
resonant cavity. From numerical solutions of the model in one dimension, we find that (i) current- 
voltage characteristics of the array have self-induced resonant steps (SIRS), (ii) at fixed disorder and 
coupling strength, the array locks into a coherent, periodic state above a critical number of active 
Josephson junctions, and (iii) when N a active junctions are synchronized on an SIRS, the energy 
emitted into the resonant cavity is quadratic with N a . All three features are in agreement with a 
recent experiment [Barbara et al, Phys. Rev. Lett. 82, 1963 (1999)]. 
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Arrays of Josephson junctions have long been studied 
both experimentally Q and theoretically Q as a poten- 
tially controllable source of microwave radiation. Most 
studies have been carried out on overdamped junction 
arrays with external loads. Typically, a dc current is in- 
jected into the array, producing ac voltage oscillations in 
each of the junctions. If all the junctions are locked to the 
same frequency, then the radiated power should vary as 
the square of the number of junctions. Overdamped junc- 
tions are usually studied, because underdamped junc- 
tions can exhibit hysteresis and chaotic behavior. How- 
ever, even overdamped arrays have proven difficult to 
synchronize: their largest experimentally achieved dc to 
ac conversion efficiency is only about 1% j|. 

Recently, Barbara et al Q achieved a 17% degree of 
power conversion in an underdamped two-dimensional ar- 
ray placed within a resonant electromagnetic cavity. In 
this case, the synchronization was achieved by an indirect 
coupling between each junction and the electromagnetic 
field of the cavity mode. The results were characterized 
by striking threshold behavior: typically no synchroniza- 
tion was achieved for arrays shorter than a certain thresh- 
old number of junctions. 

In this Rapid Communication, we present and numer- 
ically study a simple model for the dynamics of an un- 
derdamped Josephson junction array coupled to a res- 
onant cavity. This model generalizes one used recently 
to describe the energetics of such a system ||. It bears 
many resemblances to previous dynamical models, which 
either connect this array to laser action in excitable two- 
level atoms H or introduce various types of impedance 
loads to provide global coupling between junctions [[7|410[ . 
In our model, we infer the equations of motion starting 
from a more conventional Hamiltonian which describes 
Josephson junctions coupled to a vector potential |jTl| . 
Even though our model is only one-dimensional, our re- 
sults show many of the features seen experimentally j| , 
including (i) mode locking into a coherent state above a 
critical number N c of active junctions, (ii) a quadratic 
dependence of the energy on the number of active junc- 



tions above N c , and (iii) most strikingly, self- induced 
steps at voltages corresponding to multiples of the cavity 
frequency. 

We begin with the following Hamiltonian model for a 
one-dimensional array of N Josephson junctions placed 
in a resonant cavity, which we assume supports only a 
single photon mode of frequency O: 

H = H p hoton + He + Hj 

N N 

= hn{a i a+ -)+^E Cj n 2 -^E Jj cos{ lj ). (1) 
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Here, Hphoton is the energy in the cavity, He is the capac- 
itive energy, and Hj is the Josephson energy of the array. 
(V and a are photon creation and annihilation operators, 
Ecj — q 2 /{2Cj) is the approximate capacitive energy, 
and Ejj = hl c j/q is the Josephson coupling energy of 
a junction (where Cj is a capacitance, I c j a critical cur- 
rent, and q = 2|e| is the Cooper pair charge). Finally, 
-fj = ~ [(27r)/$o] Jj A • ds = (j>j - Aj is the gauge- 
invariant phase difference across a junction, where <j)j is 
the phase difference across a junction in the absence of 
the vector potential A, <I>o = hc/q is the flux quantum, 
and the line integral is taken across the junction. We as- 
sume that A arises from the electromagnetic field of the 
normal mode of the cavity. In Gaussian units, it is given 
by @,0 A(x,t) - v /(/ic 2 )/(2n) (o(t) + ot(t))E(x), 
where E(x) is the electric field of the mode, normalized 
such that J v d 3 x|E(x)| 2 = 1, V being the system volume. 
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^E(x) • ds 



(2) 



is the effective coupling to the cavity. 

The time-dependence of the operators a, a', rii, and 
<j)i follows from Eq. ([!]) together with the Heisenberg 
equations of motion ifiO = [0,H], where [..,..] is a 
commutator and O an operator. We use [a, a'] = 1, 
[a, a] = [a', a'] = z 0, and [nj,±(f>k] = ^fiSjk, with other 
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commutators equal to zero, and the relation [A, F(B)] 
[A,B]F'(B). 



We introduce the notation a = an + iai, u>, 
2u>cjOJjj, where u>cj = Ecj/h and ujjj 



pi 

Ejj/h, and 



a dimensionless natural time r = u> p t, with u> p as a 
suitable average value of u p j. For numerical conve- 
nience, we also assume that gj has the same value g 
for each junction. Then the equations of motion can 



be written (in properly scaled units) as 



— hj = 0, 

hj + (0;'^ / Qp) sin (4>j — 2a R ) — 0, a R — 05/ = 0, and 

a j + fld R — g J^ji^Jj/^p) sm {4>j ~ 25_r) = 0, where the 
dot is a derivative with respect to r. 

In order to make these exact relations amenable to nu- 
merical computation, we now replace the operators by 
c- numbers, as should be reasonable when the eigenvalues 



of 
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To introduce dissipation into the equa- 



tions of motion, we may add a term to H of the form 



ma. 2 JjY 
2 ^cc^a 



where the are random variables and Xa and p« ■* 
are canonically conjugate [Q. If the spectral density 



'h = (f ) E Q 
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-<5(w — o; Q ) = i^ay|a;|0(a; c — oj), where 
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lo c is a cutoff frequency comparable to a typical phonon 
frequency, ctj = Rg/Rj, and Rq = h/(2e) 2 , then the dis- 
sipation is ohmic |L5|] and integrating out the variables 
aia and p« leads to the usual resistively-shunted junc- 
tion equation with ohmic damping corresponding to 
a shunt resistance Rj. A driving current current can be 
included similarly by adding to H a "washboard poten- 
tial" of the form ^ X} j=i • These modifications lead 
to the following equations of motion for the 2N + 2 vari- 
ables: 



n 
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P ~ ZcCl+Aj) Q., 

a R = Qai, 

5/ = -&a R + gJ2i 



1 + Aj)sin| 



(3) 



Here, we have redefined the effective coupling as g — 
gujj/Q p , and introduced a damping coefficient Qj = 
HipRjCj, where Rj is the shunt resistance. We also in- 
troduce a disorder parameter Aj — (I C j — J C )/7 C , where 
I c is a suitable average critical current. In writing these 
equations, we have assumed that both CjRj and I c j/Cj 
are independent of j, so that each junction has the same 
damping coefficient Qj. Dissipation due to the cavity 
walls could be included similarly via a cavity Q factor. 
Note that the first two equations in (||) reduce to the 
RCSJ model in the limit of no coupling to the cavity 
(g = 0), and the last two equations to those of a har- 
monic oscillator with eigenfrequency fl. 

We have solved Eqs. (||) for the variables m, 4>i, a R and 
di numerically by implementing the adaptive Bulrisch- 
Stoer method, which is both fast and accurate We 
choose I c j, for each junction, j, randomly and indepen- 
dently from a uniform distribution between J c (l — A) and 




FIG. 1. Left scale: Current- voltage (IV) characteristics for 
an underdamped Josephson array of = 40 junctions and 
system parameters Q. — 2.2, Qj — y20, A = 0.05, and 
g = 0.001, as defined in the text. Right-hand scale: total 
photon energy in the cavity E = {a 2 R + af). Predicted volt- 
ages for the integer self-induced resonant steps (SIRS) are 
shown as dotted lines. 

I c (l + A), but for convenience assume that g is indepen- 
dent of j. We initialize the simulations with all the phases 
randomized between [0, 27r], and a R = ai = rij = 0. 
We then let the system equilibrate for a time interval 
At = 10 4 , after which we evaluate averages over a time 
interval At = 2 ■ 10 3 , using 2 16 evenly spaced sampling 
points. 

In Fig. |l|, we show the current-voltage characteristics 
for N = 40 junctions with A = 0.05 and g = 0.001, 
evaluating the time-averaged voltage from the Joseph- 
son relation, (V) = [l/Qj](YljLi 7j')- A striking feature 
of this plot is the self-induced resonant steps (SIRS), at 
which (V) remains approximately constant over a range 
of applied current. The most prominent step occurs at 
(V) / (N RI C ) = £l/Qj, but there is another, less obvi- 
ous, step at 2il/Qj. We believe the steps occur at all 
(m/n)Cl/Qj, where m and n are integers, as further dis- 
cussed below. Similar steps were seen experimentally in a 
two-dimensional array of underdamped Josephson junc- 
tions coupled to a resonant cavity As noted in Ref. 
Q, these steps are the analog of Shapiro steps in con- 
ventional Josephson junctions. They occur, we believe, 
for a similar reason: qualitatively, the variables a R oscil- 
late with frequency £1 and produce an effective ac drive 
to each Josephson junction, in addition to the dc drive 
generated by the current /. 

When we solve the system of equations (Q) numerically 
for a single junction, we find SIRS for fractions (njm) — 
1,4/3,3/2,5/3,2,5/2,3,4,.... The step width in current 
is very sensitive to g, and, indeed, we have thus far found 
the steps only for a limited range of g. For the larger 
arrays, we have not yet seen the fractional SIRS. 
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FIG. 2. Power spectrum, P(ui), of the ac voltage across the 
array, plotted versus frequency, u, at two driving currents: 
(a) and (b) I/I c = 0.492, corresponding to the first integer 
SIRS, and (c) and (d) I/I c = 0.90, slightly off a SIRS. Other 
parameters are the same as in Fig. hi (b) and (d) are the 
same as (a) and (c) except that g — 0. The vertical, dashed 
line shows the resonant frequency of the cavity. 
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Fig. fy also shows the time-averaged total energy E = 
aj) as a function of I/I c for the same array. The 
total energy in the cavity increases dramatically when 
the array is on a SIRS, and is very small otherwise. This 
sharp increase signals the onset of coherence within the 
array. 

Figure @ shows the calculated voltage power spectrum 
P(u>) = | J _ Vt t{r) exp{iu>T)d,T\ 2 , for two values of the 
driving current: I/I c = Vt/Qj [Fig. ||(a) and (b)] and 
I/I c = 0.9 [Fig. ||(c) and (d)]; all other parameters are 
the same as in Fig. El In (a), all the junctions are on the 
first SIRS and the power spectrum has peaks at the cavity 
frequency f2 and its harmonics. In (c), the array is tuned 
off the step. The power spectrum shows that the array 
is not synchronized in this case; instead, the individual 
junctions oscillate approximately at their individual reso- 
nant frequencies and their harmonics and subharmonics. 
In Fig. ||(b) and (d), we show the same case as in Fig. ||(a) 
and (c) respectively, except that the coupling constant, 
tji, is artificially set to zero. In this case, the junctions 
are, of course, independent of one another, and the re- 
sult is that of a disordered one-dimensional Josephson 
array with no coupling between the junctions. 

Next, we turn to the dependence of these properties 
on the number of active junctions, N a , in the array. The 
concept of active junction number, in the terminology of 
Rcf. |J, is meaningful only for underdamped junctions. 
As is well known, an underdamped junction is bistable 
and hysteretic in certain ranges of current, and can have 
either zero or a finite time- averaged voltage across it, 
depending on the initial conditions. In the present case, 
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FIG. 3. Left-hand scale and asterisks: Photon energy E 
in the resonant cavity when the array is current driven on 
a SIRS, plotted versus number of active junctions, N a . The 
array parameters are N = 40, Q, = 2.2, Qj = ^20, A = 0.10, 
g = 0.001, and I/I c = £l/Qj (see text). Full curve shows the 
best fit of E to the function C2N% + ciN a + co for N a > 15, the 
threshold for synchronization. Inset: E(N a ) near N a = 15, 
showing jump near synchronization threshold. Right-hand 
scale and open circles: Kuramoto order parameter, (r) (see 
text), for the same array. Dots connecting circles are guides to 
the eye. The sharp increase in (r) and the quadratic increase 
in E both start near N a ~ 15. 

N a denotes the number of junctions (out of N total) 
which have a finite time-averaged voltage drop. We can 
tune N a by suitable choosing the initial conditions, fa 
and 4>i, for each junction ||. 

We have studied the properties of the disordered array 
(A = 0.10) for N — 40 junctions, and a driving current 
I / I c = Q/Qj. This current not only lies well within 
the bistable region, but also leads to a voltage on the 
first integer SIRS. The total energy E(N a ) (normalized 
to E(6)) for this case is plotted as a function of iV a in 
Fig. ^. The active junctions are unsynchronized up to a 
threshold value N c = 15. Above this value, E increases as 
a quadratic function of N a , i. e., E = cq + c\N a + C2N?, 
where Co, c\, and ci are constants (full line in Fig. |^). 
By contrast, E is approximately independent of iV a for 
N a < N c . At N a = N c , there is a discontinuous jump 
in E by approximately a factor of 3 (see inset to figure) . 
A similar quadratic dependence above a synchronization 
threshold was also seen in Ref. 0, though for a two- 
dimensional array in an applied weak magnetic field. By 
contrast, if the system is in the bistable region, but not 
tuned to a self-induced resonant step, E does not increase 
quadratically with N a . Instead, we find E(N a ) exhibits a 
series of plateaus separated by discontinuous jumps (not 
shown in the Figure). 

To measure the degree of synchronization among the 
Josephson junctions, wc plot the Kuramoto order param- 
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eter [Q, (r) for the same parameters, as a function of 
number of active junctions, N a , (right-hand scale in Fig. 
|). (r) is defined by (r) = (|^; E^i ex P(*^)l)r, where 
(...) T denotes a time average. Note that (r) = 1 repre- 
sents perfect synchronization, while (r) = would corre- 
spond to no correlations between the different phase dif- 
ferences cj>i . As is clear from Fig. ||, there is an abrupt in- 
crease in (r) at N a = N c , indicative of a dynamical tran- 
sition from an unsynchronized to a synchronized state, 
as N a is increased past a critical value, while other pa- 
rameters are kept fixed. As expected from similar tran- 
sitions in other models ]l9p , the finite this transition is 
not inhibited by the finite disorder in the i c 's. Note that 
(r) approaches unity for large N a , representing perfect 
synchronization. This transition is the dynamic analog 
of that analyzed by an equilibrium mean-field theory in 
Rcf. S. The existence of this transition is intuitively 
reasonable, as discussed there: since each junction is ef- 
fectively coupled to every other junction via the cavity, 
the strength of the coupling increases with N a , and a 
transition to coherence should occur for sufficiently large 
N a . 

In summary, we have presented a model for a one- 
dimensional array of underdamped Josephson junctions 
coupled to a resonant cavity. We have studied the clas- 
sical limit of the Heisenberg equations of motion for this 
model, valid in the limit of large numbers of photons, 
and included damping by coupling each phase difference 
to an ohmic heat bath. In the presence of a dc current 
drive, we find numerically that (i) the array exhibits self- 
induced resonant steps (SIRS), similar to Shapiro steps in 
conventional arrays; (ii) there is a transition between an 
unsynchronized and a synchronized state as the number 
of active junctions is increased while other parameters 
are held fixed; and (iii) when the array is biased on the 
first integer SIRS, the total energy increases quadrati- 
cally with number of active junctions. All these features 
appear consistent with experiment Further study 
is underway in order to ascertain whether or not these 
features remain true of two-dimensional arrays and with 
gauge-invariant damping. 
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